function wave

%  wave equation: profiles produced using explicit method
%        diff(u,x,x)=diff(u,t,t)  for xL < x < xr, 0 < t
%  where
%        u(x,0) = f(x)  and  diff(x,t)(x,0)=-c*diff(f,x)

%  f(x) = cosine bump   

% clear all previous variables and plots
clear *
clf

% set boundary conditions and parameters
N=150;
M=272;
%M=546;
tmax=1.8;
xL=0;
xR=1;

% pick time points for plotting (by picking the index)
id(1)=1; id(2)=round((M+1)/4); id(3)=round((M+1)/2); id(4)=M+1;

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

% generate the points along the t-axis, t(1)=0 and t(M+1)=tmax
t=linspace(0,tmax,M+1);
k=t(2)-t(1);
lamda=k/h;

% calculate explicit and exact solutions
u=explicit(x,t,N+2,M+1,lamda);
ue=exact(x,t,N+2,M+1);

% plot results
% get(gcf)
%set(gcf,'Position', [1290 314 560 725]);
plotsize(560,725)
for itt=1:4

	% plot solutions at given time
	subplot(4,1,itt)
	plot(x,u(:,id(itt)),'-r')
	hold on
	plot(x,ue(:,id(itt)),'--k')
	
	% define axes and title used in plot
	ylabel('Solution','FontSize',14,'FontWeight','bold')

	% have MATLAB use certain plot options (all are optional)
	set(gca,'FontSize',14); 
	axis([xL xR -1.1 1.1]);
	box on

	% put in time value, label x-axis, and put in title
	if itt==1
		xt=0.88; yt=0.95;
		say2=['Explicit vs Exact Solution:  \lambda = ',num2str(lamda,'%3.3f'),'  (M = ',num2str(M),', N = ',num2str(N),')'];
		title(say2,'FontSize',14,'FontWeight','bold')
		text(0.06+t(id(itt)),0.7,'\rightarrow','FontSize',18,'FontWeight','bold');
	elseif itt==2
		xt=0.88; yt=0.95;
		text(0.06+t(id(itt)),0.7,'\rightarrow','FontSize',18,'FontWeight','bold');
	elseif itt==3
		xt=0.88; yt=-0.92;
		text(0.06+t(id(itt)),0.7,'\rightarrow','FontSize',18,'FontWeight','bold');
	else
		xt=0.88; yt=0.95;
		xlabel('x-axis','FontSize',14,'FontWeight','bold')
		legend(' Explicit Method',' Exact',4);
		text(1.86-t(id(itt)),-0.7,'\leftarrow','FontSize',18,'FontWeight','bold');
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
	end
	say=['t = ', num2str(t(id(itt)) ,'%3.1f')]; 
	text(xt,yt,say,'FontSize',14,'FontWeight','bold')
	
	hold off
end;

% explicit method
function U=explicit(x,t,N,M,lamda)
l2=lamda^2; k=t(2);
for i=1:N
	U(i,1)=f(x(i));
end;
U(1,2)=0; U(N,2)=0;
for i=2:N-1
	U(i,2)=U(i,1)+k*g(x(i))+0.5*l2*(U(i+1,1)-2*U(i,1)+U(i-1,1));
end;
for j=3:M
	U(1,j)=0; U(N,j)=0;
	for i=2:N-1
	  U(i,j)=l2*U(i+1,j-1)+2*(1-l2)*U(i,j-1)+l2*U(i-1,j-1)-U(i,j-2);
	end;
end;

% exact solution
function U=exact(x,t,N,M)
a=0.09;
tt=2-a;
for j=1:M
	for i=1:N
		U(i,j)=f(x(i)-t(j))-f(x(i)+t(j)-tt);
	end;
end;

% subfunction f(x)
function q=f(x)
x0=0.09;
q=0;
if (x>0)&(x0)&(x